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There is increasing interest in many-body perturbation theory as a practical tool for the calculation 
of ground-state properties. As a consequence, unambiguous sum rules such as the conservation 
of particle number under the influence of the Coulomb interaction have acquired an importance 
that did not exist for calculations of excited-state properties. In this paper we obtain a rigorous, 
simple relation whose fulfilment guarantees particle-number conservation in a given diagrammatic 
self-energy approximation. Hedin's Go Wo approximation does not satisfy this relation and hence 
violates the particle-number sum rule. Very precise calculations for the homogeneous electron gas 
and a model inhomogeneous electron system allow the extent of the nonconservation to be estimated. 
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I. INTRODUCTION 

Many-body perturbation theory is a powerful method 
for studying interacting electron systems, because the 
partial summation of self-energy diagrams allows an ef- 
ficient and systematically conveijffing description of the 
dominant scattering mechamsmsu In solid-state physics, 
Hedin's GW approximatiorfl includes dynamic screening 
in the random-phase approximation and has been apt. 
plied with great success to a large range of materials.El 
While calculations have long focustd on electronic ex- 
citations, such as band structurespQ that are not nor- 
mally accessible by variational mean-field schemes, there 
is now increasing interest in using many-body perturba- 
tion theory also ta obtain ground-state, ptpperties like 
the charge densityQ or the total energylZrt3 in order to 
circumvent well-known limitations o£^tandard approxi- 
mations in density-functional theory.llil Unlike the calcu- 
lation of excited states, which are given immediately by 
the pole structure of the spectral function, this generally 
requires a multidimensional integration over the hole part 
of the Green function, as Hi|Galitskii and Migdal's expres- 
sion for the total energyJ13 As a consequence, sum rules 
that could hitherto be ignored have gained new promi- 
nence. The most important of these is the conservation 
of particle number, i.e., the requirement that the integral 



N = 



(1) 



over the diagonal elements of the Green function G equals 
the true number of electrons. Here a denotes the spin 
variable and 77 is a positive infinitesimal that forces the 
frequency contour to be closed across the upper complex 
half plane. i— . 

In a seminal paper Baym and KadanofiliS investigated 
the evolution of nonequilibrium Green functions and 
derived a set of symmetry relations for diagrammatic 
many-body approximations that guarantee the overall 



conservation of particle number, total energy, and mo- 
mentmn under time-dependent external perturbations. 
Baymliil later showed that a self-energy satisfying all 
of these relations can be represented as the derivative 
S = S^/SG of a generating functional <&, and that in 
this case the Green fimction obtained self-consistently 
from Dyson's equationlis moreover yields the exact par- 
ticle number. In particular, this applies to the fully self- 
consistent GW approximation, in which both the Green 
function G and the screened Coulomb interaction W are 
dressed by self-energy insertions in accordpace with a 
self-consistent solution of Dyson's equation.Ej However, 
it has since become clear that $ derivability is a suffi- 
cient but not a necessary requirement for the fulfilment 
of the particle-number sum rule. For instance, the par- 
tially self-consistent GWq approximation, in which only 
the Green function is updated self-consistently but the 
screened Coulomb interaction remains undressed, is not 
$ derivable but nevertheless produces the exact particle 
number.ta On the other hand, without self-consistency 
even in the Green functicm, the particle number is not, 
in general, given correctlyJ13 but this computationally ef- 
ficient GqWq approach still remains the preferred method 
for most practical applications. 

More complicated self-energy exfpessions like the cu- 
mulant expansiontSI or the T matrixll3 have already been 
successfully applied to solids, leading to an improved de- 
scription of satellite resonances. Like the GW approxi- 
mation, these schemes are typically implemented without 
full self-consistency and are hence not $ derivable. In or- 
der to avoid expensive numerical tests in such situations, 
it would be desirable to have clear diagrammatic criteria 
for the fulfilment of the particle-number sum rule that 
could be checked a priori without actual calculations. 

Unfortunately, Baym's proof, which is based on Lut- 
tinger's examination of the exact theorycEl and deter- 
mines the volume of the Fermi sea directly, cannot eas- 
ily be extended, because it relies explicitly on the exis- 
tence of the generating functional <&. We therefore take 
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a different approach by describing the switching on of 
the Coulomb potential as a time-dependent process that 
connects the noninteracting and the corresponding inter- 
acting electron system on a finite time scale. In this way 
we can examine the differential conservation laws and 
deduce a diagrammatic symmetry relation for particle- 
number conservation before taking the adiabatic limit. 
The theoretical framework is developed in Sec. The 
non-self-consistent GqWq approximation, which violates 
this symmetry relation and does not conserve the parti- 
cle number when the interaction is switched on, deserves 
special attention owing to its pre-eminent role in practi- 
cal implementations. In Sec. Ill we therefore present very 



precise numerical calculations of the particle number for 
the homogeneous electron gas and a model inhomoge- 
neous system in order to assess the quantitative devia- 
tion. Finally, in Sec. IV we summarize our conclusions. 
Atomic units are used throughout. 



II. PARTICLE-NUMBER CONSERVATION 

In order to connect the interacting electron system 
with the corresponding noninteracting system, whose 
properties are supposedly known exactly, we consider the 
Hamiltonian 



H{t) = Ho 



(2) 



The one-body part Hq contains the kinetic energy as well 
as the external potential V^xt, and the Coulomb interac- 
tion Hi is switched on exponentially with e > 0. At large 
times, both in the past and in the future, the Hamiltonian 
reduces to Hq, which constitutes a solvable problem. The 
noninteracting Green function Gq is readily constructed 
from the solutions of the single-particle Schrodinger equa- 
tion and yields the correct number of particles. On the 
other hand, at t = the full Coulomb interaction is ef- 
fective, and the Green function is defined aM 



G{x,x') 



{^T[i:{x)i;\x'm) 



(3) 



where the shorthand notation x = (r, a, t) indicates a set 
of spatial, spin, and temporal coordinates, \^) denotes 
the ground-state wave function of the interacting elec- 
tron system in the Heisenberg picture, and T is Wick's 
time-ordering operator that rearranges the subsequent 
symbols in ascending order from right to left with a sign 
change for every pair commutation. Furthermore, ?/'^(a;') 
and ■0(2^) represent the electron creation and annihilation 
operator in the Heisenberg picture, respectively. 

The unknown many-body wave function |^) evolves 
from the noninteracting ground state |^'o) and can for- 
mally be expressed as j^*) = ?7e(0, — oo)|5'o), whereH 



i/=0 



dU 



-e(|tl| + ---+|t,.|) 



represents the time-development operator. With this def- 
inition the Green function may be rewritten as 



G{x, x) 



{^o\T[S,i;{x)tlj\x')]\^o) 

(*0|^e|*0) 



(5) 



with S'g = U^UiD^ — oo). At this stage the Gell-Mann and 
Low theoremlEl asserts that it is, in general, permissi- 
ble to take the adiabatic limit e ^ 0. However, in the 
following we continue to perform a time-dependent per- 
turbation analysis for finite e and only take the adiabatic 
limit after establishing the conservation criteria that ap- 
ply during the transition. The time-ordered products 
in Eq. (|^) may be evaluated in the usual way by in- 
voking Wick's theorem,E3 because the exponentials are 
scalar functions and commute with the field operators. 
Hence the perturbative treatment generates the standard 
series of cojinected and topologically distinct Feynman 
diagramsJHj made up of the noninteracting Green func- 
tion Go and the two-body Coulomb potential, but the 
latter now acquires an additional prefactor and is given 

by 



v{x, x') 



r — r' 



(6) 



The formal identity of the perturbation expansion in the 
time-dependent and the adiabatic, time-independent case 
is a crucial result that forms the basis of our discussion 
in this section. 

For an analysis of the conservation properties we now 
follow Ref. O and write the perturbation series as 



Gq "^{x, xi)G{xi, x') dxi 

5{x — x') — i v{x,xi)G2{x,xi;x' ,x^)dxi , (7) 



invoking the two-particle Green function G2 . The super- 
script a;"*" indicates that a positive infinitesimal is added 
to the time variable to ensure the proper ordering. An 
equivalent form is the adjoint equation of motion 



G{x, xi)Gq^ {xi , x') dxi 

5{x'~x)~i G2{x, xi; x' , Xi)v{x' , Xi) dxi . (8) 



X 



T[Hi{h) ■ ■ ■ Hi{t,)] 



(4) 



The inverse noninteracting Green function is identical to 
the operator 

= (^-i^ + \v' -V,^,{v')y{x-x'), (9) 
so that after substracting Eq. (||) from Eq. (Q) we obtain 
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M I: + T?:^ + ^ (V + V) • (V - V) 



G{x,x') 



ydt dt' J ' 2 
[Kxt(r)-Fe,t(r')]G(x,a;') (10) 

—i I [v{x, xi) — v{x , xi)] G2{x, x\\x , ) dx\ . 



When we set x' — x^, the terms on the right-hand side 
cancel, while the left-hand side reduces to the differential 
conservation law for the particle number 



dn(y, t) 
dt 



V-j(r,t) = 



(11) 



with the electron density n{r,t) = —i^^G{x,x^) and 
current j(r,i) = -i E.KV - V')G(a:, x')],'=.+ . Thus 
whenever Eqs. and (||) are satisfied simultaneously, 
the total particle number is conserved while the interac- 
tion is switched on. This does not depend on the value of 
e and, in particular, remains true in the adiabatic limit, 
which can now be taken, turning G into the equilibrium 
Green function of the interacting electron system. 

By multiplying Eqs. and (^) with G from the left 
and right, respectively, and then subtracting one from 
the other, their mutual consistency can be stated in the 
more convenient form. 



G{x, X2)v{x2,Xi)G2ix2,Xi;X , Xi ) dxi dX2 

G2ix,xi;x2,Xi)v{x2,xi)G{x2,x') dxidx2 , (12) 



that may easily be verified by visual inspection of a given 
diagrammatic approximation for the two-particle Green 
function. Evidently it is, the same criterion as derived 
by Baym and KadanofitJ for particle-number conserva- 
tion under time-dependent external perturbations. This 
is no coincidence, of course, because the termwise can- 
cellation of diagrams on the right-hand side of Eq. ( |l0| ) 
is of purely topological origin and does not depend on 
the mathematical properties of the constituent propaga- 
tors. Hence it is inconsequential whether, as in Refs. |l3| 
and ^4|, Go contains a time-dependent perturbation while 
the interaction is constant or, as in the physical situation 
considered here, the noninteracting Green function is in- 
variant under temporal translations while the Coulomb 
potential instead acquires a time-dependent prefactor. 

As an example we now consider the GW approxi- 
mation. The self-energy, when applied with full self- 
consistency, is given by 



S(a;,x') = iG{x,x')W{x+,x') 



(13) 



where the screened Coulomb interaction W takes the 
mathematical form of the random-phase approximation 
but is evaluated using the dressed Green function self- 
consistently derived from Dyson's equation. The dia- 
grammatic representation of S is shown in Fig. |l|(a) . The 
corresponding two-particle Green function is obtained by 
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FIG. 1. Diagrammatic representation of the self-energy 
T.{xi,X2) and the corresponding two-particle Green function 
G2{xi, xs; X2, X4) in (a) the fully self-consistent GW approxi- 
mation, (b) the partially self-consistent GWo approximation, 
and (c) the non-self-consistent GqWq approximation. 



comparing Dyson's equation with the equation of motion 
(0), which yields the identity 



~i J v{x,xi)G2{x,xi]x' ,x^) dxi 

= Vh{x)G{x, x') + T,{x, xi )G{xi , x') dxi 



(14) 



where Vh{x) — Jv{x,xi)G{xi,x^) dxi indicates the 
Hartree potential. The two-particle Green function cor- 
responding to the GW approximation for the self-energy 
is also displayed in Fig. |l|(a). It is easily seen that it sat- 
isfies the symmetry relation (^2|), which is essentially a 
horizontal left-right symmetry for the building blocks of 
G2, and hence conserves the total particle number when 
the Coulomb interaction is switched on. Of course, this 
result also follows from the existence of the generating 
functional <i>.Efl 

The partially self-consistent GWq approximation 



'S{x, x) — iG{x, x')Wo{x'^ , x') , 



(15) 



in which the screened Coulomb interaction Wq is evalu- 
ated with the noninteracting Green function Gq, is not 
<i> derivable, which would require an additional vertical 
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mirror symmetry 6*2(3:1, X3; X2, 0:4) = 6*2(2:3, cci; 0:4, 0:2) 
in the diagrammatic structure of the two-particle Green 
function that has been lost in the transition from full 
to partial self-consistency. Nevertheless, G2, shown in 
Fig. ^b), still obeys the consistency relation ( p^ and 
hence guarantees the correct total particle number, as 
previously ccmfirmed by explicit integration of the spec- 
tral function.tj In contrast, the non-self-consistent GqWq 
approximation 



Our concern is the evaluation of the particle-number 
difference. 



T.{x,x') = iGo{x,x')Woix+,x') 



(16) 



leads to a two-particle Green function with lower internal 
symmetry, displayed in Fig. 0(c), that no longer satisfies 
Eq. (p^, implying an incorrect total particle number. 
The quantitative deviation is investigated in the follow- 
ing section. In a similar manner, the conservation prop- 
erties of other diagrammatic self-energy approximations 
are easily established by an inspection of the underlying 
two-particle Green function. 



III. NUMERICAL RESULTS 

In the previous section we proved that the GW and 
GWq approximations conserve the particle number for an 
arbitrary electron system when the Coulomb interaction 
is switched on, in contrast to GqWo- This, coupled with 
their superinrj^jerformance in ground-state total-energy 
calculations ,LrEj might be thought to suggest that the 
GqWq approach is useless if one is interested in ground- 
state properties. However, a many-body calculation at 
only the GqWq level is already sufficient to correct typ- 
ical limitations of mean-field density-functional theories, 
such as their inaccuracy in highly inhomogeneous-Sjystems 
or their failure to describe van der Waals forces .cJ More- 
over, the Green function arising from a GqWq calculation 
may be used as |iiiput in the variational Luttinger and 
Ward functional ,E3 and prospective calculations suggest 
that this_js an excellent approach for calculating total 



energiesJij Since these methods are more amenable to ap- 
plications in complex systems than the fully or partially 
self-consistent GW approximations, it is important to de- 
termine whether the underlying violation of the particle- 
number sum rule in the GqWq framework is small enough 
to be safely ignored. 

There are some indications that such an error is in- 
deed fairly small Jpr the homogeneous election gas at 
metalHc densities a Hubbard model system,tZl and typ- 
ical semiconductors .□ Here, bearing in mind that many- 
body total-energy calculations are intended to be used in 
extreme situations where standard implementations of 
density-functional theory fail, we present numerical re- 
sults for thin jellium slabs, whose most relevant feature 
is the strong inhomogeneity of the electron-density pro- 
file, as well as for the homogeneous electron gas over a 
wide range of densities. 



6N = — 

TT 



dujiY [G[u) ~ Gq{lo)] 



(17) 



where tr denotes the spatial trace (we omit the explicit 
spatial variables for clarity and also consider only spin- 
unpolarized systems). Since both G and Go behave as 
I/ld for large frequencies, we can apply Cauchy's theorem 
and write Eq. ( p7| ) alternatively as 

1 /■+°° 

SN^- dLutr[G{^l + iuJ)^GQ{^lQ+iuJ)] , (18) 

J-00 

where /i and fiQ are the chemical potentials of the inter- 
acting and the noninteracting system, respectively, which 
correspond, by definition, to the position of the pole of 
the Green function at the Fermi surface. As the charac- 
teristic sharp structure of G{oj) (quasiparticle peaks and 
satellites) does not appear in the analytic continuation 
G(/i + iuj), Eq. ( ^ is preferred for numerical integra- 
tion. We hepce follow some of the ideas suggested by 
Rojas et alzH and work exclusively in an imaginary time 
and frequency representation. An accurate evaluation of 
Eq. ( p^ ) furthermore requires a treatment of the high- 
frequency tails of G, which can be done easily with the 
numerical procedures described in Ref. 
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For the homogeneous electron gas, an analytic expres- 
sion exists for the noninteracting Green function Go (r, zt) 
in real space and imaginary time,c3 while the screened 
Coulomb interaction Wo{k,iuj) in the random-phase ap- 
proximation is given analyticalUi,in reciprocal space by 
the dynamic Lindhard function£3 The evaluation of the 
self-energy according to S(r, ir) = iGo{r,iT)WQ{r,iT) 
therefore only requires the numerical Fourier transform 
Wo{k,iuj) WQ{r,iT). It is this largely analytic ap- 
proach that makes the present calculation especially pre- 
cise. 

At this stage we remark that the self-energy given by 
Eq. ( p6| ) has the same analytic structure as the under- 
lying Green function Gq, i.e., the poles of S(w) are lo- 
cated in the upper (lower) complex half-plane for energies 
smaller (larger) than ^0 = ^kp, where fcp denotes the 
Fermi wave vector. As a consequence, an inconsistency 
arises because the true self-energy should have a polar 
structure identical to the interacting Green function with 
the chemical potential fi. The self-energy must therefore 
be appropriately shifted along the real frequency axis. In 
the imaginary time/frequency representation, this shift is 
automatically included in the backward transform 



I](/i + iuj) = 



+00 



dr E(iT)e 



(19) 



The calculation of S(/i -I- iuj), therefore, does not require 
an advance knowledge of /i, which can now be obtained 
from the relation fi = fiQ + T,{k-p, /i). 
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FIG. 2. Violation of the particle-number sum rule for the 
homogeneous electron gas in the Go Wo approximation. The 
relative error in the density is always negative and of the order 
of 0.1% in the range of metallic densities. 



Finally, the interacting Green function is calculated in 
reciprocal space according to 



G{k, /i + iijj) 



1 



(20) 



In the same representation, the noninteracting Green 
function is given by 



1 



(21) 



so that the density variation is readily obtained from 



6n 



(2^ 



+ 00 



duo 



[G{k, fi + iLu) — Go{k, Ho + i^)] 



(22) 



In Fig. ^ the relative devation Sn/no from the exact den- 
sity is displayed as a function of the Wigner-Seitz radius 
Tg. In the high-density region < 1.8 the particle num- 
ber is slightly overestimated (< 0.01%), while it is un- 
derestimated for lower densities. In the range of metallic 
densities this underestimation is of the order of 0.1%, 
but the error becomes increasingly important in the di- 
lute limit (-1.7% for = 10 and -6.1% for = 20). 

As pointed out above, it is also of interest to investi- 
gate the error resulting from the GqWq method in the 
total number of particles for a strongly inhomogeneous 
system. The model we have chosen is a thin jellium slab 
with a background density uq — {^TTr^)~^ and width 
L. The slab is bounded by two infinite planar walls, 
so that, if charge neutrality is assumed, the system is 
fully characterized by the lengths and L. In this case. 
Go corresponds to the Kohn-Sham system obtained self- 
consistently with the local-density approximation (LDA) 
for the exchange-correlation potential T^c, as is typically 
done in practical ab initio calculations. 



With z chosen as the coordinate perpendicular to the 
planar walls, the translational symmetry of the system 
in the xy plane allows an efficient semianalytic evalua- 
tion of the relevant propagators. The screened Coulomb 
interaction is given by Wo — ^q^v, where eo denotes 
the dielectric function in the random-phase approxima- 
tion. The latter is calculated as eQ{k^iuj)af3 in the ba- 
sis Ca(2)exp(ik • p)/\/S. Here Ca{z) is a set of co- 
sine functions, k = {kx,ky) and p — {x,y) denote 
the two-dimensional momentum and the position vec- 
tor in the xy plane, respectively, and S is the slab sur- 
face. The matrix elements can be calculated anah^tically 
in terms of the scalar products (CQ0n|0m),Oi^ where 
^„(2:)exp(ik-p)/-\/5 are the single-particle eigenstates of 
the Kohn-Sham Hamiltonian /iks- The matrix elements 
v{k)ai3 of the Coulomb potential are likewise obtained 
analytically. The screened Coulomb interaction is then 
easily calculated by a matrix inversion for each value of fc, 
and the real-space representation is given by expanding 



/(Pk f 



2^ 



xCa(z)C/3(^')W'o(fc,Zw)a/3 . (23) 

The Green function Go(p, z, z';iT) is readily calculated 
from the Kohn-Sham eigenstates, and by employing 
Eq. ([l^ ) we obtain the self-energy in real space and 
imaginary time as well as, eventually, its representation 
S(A:, H+iuj)nm in the Kohn-Sham basis set. The presence 
of infinite confining walls implies a quick convergence 
with respect to the number of cosine and Kohn-Sham 
wave functions used in the calculation. The convergence 
is further accelerated by the analytic treatment of the 
asymptotic time and frequency tails of all operators. 

The Green function is calculated in the basis of Kohn- 
Sham eigenstates according to 

G{k, ji + iuj) = [iuj — ft,KS (fc) ^ S(fc, + iuj) 

+ V^cik)+tir^ (24) 

by a matrix inversion in the indices nm. Finally, the 
variation of the number of particles per surface unit is 
given by 



SN 



E 



d'k 

(2^)2 



dw 



[G{k, /i + ius) 

r. 



r 



(25) 



where we have used the invariance of the trace with re- 
spect to any wave-function representation. 

In Fig. H we plot the relative deviation of the parti- 
cle number SN/N in the GqWq approximation for sev- 
eral configurations of the model system, keeping the ex- 
act number of particles per surface unit n^^ = uqL = 
L/(|7rrg ) constant. The limit L ^ thus corresponds to 
a two-dimensional (2D) homogeneous electron gas with 
density n^^. Over the wide variation of the degree of 
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FIG. 3. Relative violation of particle number in the GoWo 
approximation for thin jellium slabs of fixed 2D density 
= 3/iiT as a function of their thickness L (and the cor- 
responding 3D density parameter rs). A typical error bar is 
reported. 



homogeneity shown in the figure, it is seen that SN/N 
remains of similar magnitude as in the homogeneous case 
(^ 0.2%). This observation remains true for other 2D 
densities inside the range [0.1,1]. 



IV. CONCLUSIONS 

In this paper we have rigorously obtained a general 
criterion which allows, by simple inspection, to verify 
whether a diagrammatic self-energy approximation satis- 
fies the particle-number sum rule for an interacting elec- 
tron system. As an application, we have demonstrated 
that the so-called GqWo method does not yield the cor- 
rect particle number, generalizing the conclusions of a 
previous analytic study for a Hubbardanodel Hamilto- 
nian defined only on a discrete lattice.llZl Thus this lim- 
itation of the Go Wo approximation has been fully con- 
firmed for arbitrary electron systems. By performing a 
very precise integration of the spectral function, we have 
furthermore calculated the size of the error in the GqWq 
particle number in two simple, but very distinct, fami- 
lies of electron systems. The error becomes large only 
outside the range of densities of physical interest. 
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